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Abstract 

A new Monte Carlo approach is proposed to investigate the fluid-solid phase transition of the 
polydisperse system. By using the extended ensemble, a reversible path was constructed to link 
the monodisperse and corresponding polydisperse system. Once the fluid-solid coexistence point of 
the monodisperse system is known, the fluid-solid coexistence point of the polydisperse system can 
be obtained from the simulation. The validity of the method is checked by the simulation of the 



'^ ' fluid-solid phase transition of a size-polydisperse hard sphere colloid. The results are in agreement 

O , with the previous studies. 
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Melting and freezing transitions are among the most important phenomena in condensed 
matter physics. The determination of the equihbrium phase diagram of a specific material is 
the key point in the understanding of the equilibrium properties of the material. The hard 
sphere system(HSS), which is a representative model of a large class of colloidal suspensions, 
has been extensively investigated. It is well known that HSS undergoes an entropy driven 
first-order phase transition from a disordered fluid to an ordered solid as the packing fraction 
increases |l|. However, in a true colloidal system the particles are inevitably polydisperse in 
particle sizes, which can be measured by the ratio of the standard deviation to the mean of 
the diameter, S = ^ _ "^ . It has a remarkable effect on the thermodynamic and dynamic 
behaviors of the HSS |2|, yl, U, la, la, [Zl]. The influence of the polydispersity on the fluid- 
solid phase transition of the HSS has been studied experimentally [2], and theoretically by 
computer simulation |3| , density functional theory Sl, l9| , moment free energy method [4| and 
other simplified theorie s [id . 111. Il2l|. Experimentally, crystallization process depends on the 



nonequilibrium effects [l3l, ll^l , as a result the observed phase behaviors deviate more or less 
from the equilibrium case. However, the theoretically obtained phase diagram is usually the 
equilibrium one. Therefore, the simulation study of freezing transition of the polydisperse 
HSS is essential to check the validity of the theory and to deeply understand the equilibrium 
phase behaviors of the system. 

When studying the fluid-solid coexistence of the polydisperse HSS by computer simu- 
lations, there exist two basic difficulties which also exist in the monodisperse case: one is 
that the coexisting solid forbids the particle insertion; and the other is that the back and 
forth tunneling event between the liquid and solid phases never happens within the sim- 
ulation time scale. Therefore, some usual methods, by which the phase transition of the 



polydisperse soft spheres 15l, ll6| and the nematic-isotropic transition of polydisperse liquid 
crystal 17| can be investigated effectively, are not suitable for the problem. To the best of 
our knowledge, Gibbs-Duhem integration algrithm[3| is the only available assumption- free 
method to determine the phase behavior of polydisperse HSS, where the coexistence curve is 
obtained by integrating the Clausius-Clapeyron equation from the monodisperse hard sphere 
coexistence states. However, the method is not very robust and lacks built-in diagnostics, 
and no other simulation data are available for comparison. Thus it is desirable to develop 
such an approach by which the phase transition can be determined from the rigorous free 
energy calculation. One possible scheme is to directly estimate the free energy difference 



between the stable polydisperse hard 



sn 



lere fluid and solid phases. In principle, the semi- 



grand ensemble lattice switch method 18|, ll9| is a feasible one. However, the implementation 
of a simulation along this line is nontrivial. First, the free energy difference between the 
fluid and the solid phases is huge, second, the system need to traverse a enormous entropy 
barrier ("gate" state), and third, the symmetry of fluid phase is not the same as solid phase. 



The other possible scheme is the so called referenced-state method [20|, l2l| . For the polydis- 
perse HSS the best reference states are just the coexisting fluid and solid of monodisperse 
HSS. Unfortunately, in the simulation it is often difficult to find a reversible path linking the 
monodisperse reference state to the polydisperse system under consideration. The reason is 
that the ensemble for simulating the polydisperse system(semigrand ensemble) is different 
from that for the monodisperse system(canonical ensemble). 

In the present work, we use the concept of extended ensemble to construct a reversible 
path connecting the polydisperse target state to the monodisperse reference state, along 
which the free energy difference between polydisperse and monodisperse system can be esti- 
mated easily. By the method we reinvestigate the fluid-solid phase transition of polydisperse 
HSS. 

The semigrand canonical ensemble (SCE) is the best frame to investigate the phase be- 
havior of a polydisperse HSS. In this ensemble the independent variable is the chemical 
potential difference function, the total number of particles is flxed and the number of par- 
ticles of each size is permitted to fluctuate. A continuous composition distribution can be 
realized on average. Generally the composition distribution can not be prescribed in ad- 
vance unless we know the conjugated chemical potential difference function. Fortunately, 
the chemical potential difference function can easily be obtaine d by the recently developed 



SNEPR method, which was described in detail in reference [18|, |27|] . In all simulation stud- 
ies of polydisperse systems the polydisperse parameter(diameter of particles here) is always 
discretized. For the purpose of describing our method more easily, we present the semi- 
grand canonical ensemble in the case of the discrete diameters. The diameter of particles 
is assumed to take M + 1 discrete values (Tj(A) = am + {i — y)'M' ^^^^ ^ = 0, 1, ■ ■ ■ , M, 
A = amax - cTmin, cTm = {(^max + crmm)/2, (Tmax and amin are the maximum and minimum 
of the allowed diameters of the particles, respectively. In the limit of monodisperse case, 
(7max = o'min = o'm, A = 0, and all (jj's are the same. However, we can still regard the 
particles with different i's as belonging to different species. We will see that this concept is 



useful for later construction of reversible paths between a monodisperse and a polydisperse 
system. 

By introducing the excess chemical potential relative to ideal gas fiexi^'i) = fJ'{o'i) — 
A;Tln( — y-), the partition function T for the semigrand canonical ensemble is written as 

^ " mX^^ia ) 5Z ■ ■ ■ 5^ ^^ ^ ^^P \|3^{^^ex{d^) - /iex(ffr)) \ ■ (1) 

Here, (3 = 1/kT, a,, is the diameter of an arbitrary chosen referenced component, A(o"r) = 
h/{27rmrkTy^'^ is the thermal wavelength of the referenced component, da is the diameter 
of the ath particle, and Zn is the canonical configuration integral 

Zn= f ■■■ f e-^^n^r^. (2) 

Jri Jtm a=l 

The partition function T is related to the semigrand canonical free energy Y through 

r = -A;TlnT(iV,V^,T,{A/ie.(a,)}), (3) 

here A/iex(o-j) = Atex(o-j) - /iex(o-r)- 

In order to establish a reversible path that links the polydisperse system to the monodis- 
perse referenced system, we employ an extended semigrand canonical ensemble in the simu- 
lation. In the extended ensemble the range of particle sizes A = (Jmax — o'min is regarded as 
a new ensemble variable. The extended ensemble is composed of the semigrand canonical 
ensembles with different range of particle sizes A. Each value of A specifies a macroscopic 
state of the extended ensemble. In what follows, it is convenient to set the diameter of the 
referenced species a,. = o"m,- With the additional ensemble variable A, the partition function 
of the extended semigrand canonical ensemble is defined as 

r(iV, V, {A/Xex(^.)}) = Yl ^(^' ^' {A/^ex(^.)}, A), (4) 

{A} 

here A takes a series of discrete values ranged from A = to A = Amax, corresponding 
to the monodisperse reference system and the polydisperse target system, respectively. The 
{Afiexio'i)} are kept the same for all values of A. For A = 0, the monodisperse case, 
these {Afiex{<^i)} give a distribution of particles among different species with the same 
diameters, while for A = A^ax the {Afj,ex{cri)} give the required distribution of particles 
of different sizes. Thus, the monodisperse system can be transformed into the polydisperse 



system by the variable A. For the sake of convenience, the mid-point of the particle sizes, 
o"m = {o-max + o'min)/'^, IS fixed when changing A. In the simulation the A moves are 
performed by resizing all particles, however, the partition of the particles among species 
remains unchanged even though the particle sizes are changed. The semigrand canonical 
free energy difference between the macroscopic states A = A^ax and A = is written as 

Pr(0) 



AF = ln 



(5) 



where Pr(A) is the probability that the system is in the macroscopic state A. In the 
extended semigrand canonical ensemble the probability can be obtained from simulation by 



the flat histogram methods 23, l25|, |26|]. In the simulation the trial moves include particle 
translations, resizing particle operations(breathing move) and changing A operations. 
The semigrand canonical partition function with A = can be written as 

T(iV,V^,{A/z,.(a,)},0) = ^^C(iV,T,{A/ie.(^.)}) (6) 

m 

Where the configuration integral Z^ is independent of the composition distribution when 
A = 0. And C{N, T, {Afiex{o'i)}) is a summation in the particle species, independent of the 
volume. Equations (jS]) and (jnD bridge between the semigrand canonical partition function 
of polydisperse system and the canonical partition function of monodisperse system. 

As was well known, the fluid-solid phase transition of monodisperse HSS can be deter- 



mined by looking for the equal- weight double-peak structure in the volume histogram 22 1 . 
In actual computation we use a more tractable criteria "equal peak height" to locate the 
transition point [23], which differs from equal-weight in the finite-size effect. They give the 
same result in the thermodynamic limit. If the number of particles of the coexisting fluid 
is equal to that of the coexisting solid, according to the "equal peak height" condition their 
canonical partition functions will satisfy the following relation 

Z^exp(-/3P^V7) = Z:exp{-f3P'^V:), (7) 

here, P'^ is the coexistence pressure, V? and V/ are, respectively, the volume of the monodis- 
perse coexisting fluid and solid, and ZJ and Z^ are the canonical partition functions of 
monodisperse coexisting fluid and solid, respectively. Substituting ([6]) in ([7]), we get 

TjiV;, {AfxU^,)}, 0)exp(-/3P^V7) = T,(K^ {A/i,.(a,)}, 0)exp(-/3P^K'), (8) 
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FIG. 1: The flow diagram of the simulation of polydisperse fluid-solid coexistence. The acronym 
in the small boxes correspond to: MCF and MCS: Monodisperse Coexisting Fluid and Solid; PFO 
and PSO: Polydisperse Fluid and Solid with A = 0; PTF and PTS: Polydisperse Target Fluid and 
Solid with A = Amax] PCF and PCS: Polydisperse target Coexisting Fluid and Solid. The dotted 
lines describe the coexistence of the fluid and solid phases. 



where C{N,T,{Afiex{<^i)}) is eliminated for it independent of the state of system. Thus, 
in the polydisperse system with A = the fluid phase is connected to the solid phase. For 
the polydisperse system with A = A^ax that we consider, equations ([5]) and ([8]) establish a 
relation between the fluid phase and the solid phase. 

Therefore, once the phase transition point of the monodisperse system is located, com- 
bining equations ([5]) and ([8]) the semigrand canonical free energy difference between the solid 
and fluid phases of the polydisperse target system with given Amax and A/iex(o"j) can be 
calculated by the extended semigrand canonical ensemble simulation where the A is taken 
as the ensemble variable. In real calculations the monodisperse coexistence point was taken 
from literature where the fluid volume fraction is 0.492 and the solid volume fraction is 
0.543, respectively. The remaining task to locate the transition point of the polydisperse 
target system is trivial. The polydisperse target fluid and solid phases are treated as two 
new referenced states, then the distribution of the volume of the polydisperse target fluid 
and solid phases can be obtained separately by performing another extended semigrand 
canonical ensemble simulation, taking the volume V as the ensemble variable. From such 
extended ensemble simulations the T{V) can be calculated directly by the flat histogram al- 
gorithm. When the volume distribution T{V)exp{—f3PV) shows two peaks of equal height, 
the fluid-solid phases coexistence occurs. The flow diagram of the calculation procedure is 
shown in figure [H 

To test the validity of our method, we use it to locate the fluid-solid transition point 
of the polydisperse hard spheres. Our system is composed of 256 size polydisperse hard 
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FIG. 2: Left: the chemical potential difference as a function of the particle diameter. Right: 
the corresponding composition distributions of the coexistence fluid phase(open circles) and solid 
phase(filled circles). The dotted lines denote the prescribed composition distribution of the initial 
target polydisperse crystal, here the upper one with a size polydispersity of 0.0163, the middle one 
with a size polydispersity of 0.0252, and the lower one with a size polydispersity of 0.0507. 



spheres contained in a periodic box. The crystal structure under consideration is face- 
center-cubic (fee), since our previous calculation [l8| showed that fee phase is still the most 
stable structure of polydisperse hard sphere crystal. The composition distribution of the 
initial target polydisperse solid linked to the referenced monodisperse solid is the truncated 
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FIG. 3: The probability distribution of the dimensionless volume of the system at coexistence 
pressure, the imposed chemical potential difference function is the upper one in Fig [2j The di- 
mensionless volume is defined as V* = V/N{a)'^, here a is the average diameter of particles of the 
initial polydisperse target crystal. 

Schultz function. The conjugated chemical potential difference function is evaluated by the 
SNEPR method. The solved chemical potential difference function is saved and then used 
to perform the rest of the simulation. The initial target fluid is produced by using the solved 
chemical potential difference function. Of course, it is possible to prescribe the composition 
distribution of the initial target fluid rather than the initial target solid. Previous studies 
by Bolhuis et al used the quadratic form for the chemical potential difference function |3|, 
the corresponding composition distribution is not the same as ours. However, for the same 
polydispersity 6 we expect that our simulation will give qualitatively similar results. 

Using the SNEPR method we calculated the chemical potential difference functions for 
three different truncated Schultz distributions in the fee solid phases. Each distribution 
has different polydispersity, and the maximum one is about 6 ~ 5% because at higher S 



the crystal may be unstable 
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]. Figl2] shows the chemical potential difference 
functions and the relevant composition distributions of the coexisting phases. We see that 
the particle size distribution of the coexisting fluid phase significantly differs from that of the 
coexistence solid phase, this is the phenomenon named as fractionation, the fractionation 
effect increases with the polydispersity. Comparing with the coexisting solid, the fluid 
has larger polydispersity and lower volume fraction. These results are consistent with the 



previous results 



291 ]. We also noted that the particle sizes distribution of the coexisting 
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FIG. 4: The phase diagram of the polydisperse hard sphere system in the {5,p) plane. The dotted 
Unes denote the boundaries of the fluid and solid phases. The solid lines denote the tie lines 
connecting the two coexisting phases. 



solid differs from that of the prescribed target solid, this is because the volume of initial 
target solid is not equal to that of the coexisting solid, or we can determine directly the 
cloud points(the coexisting phase with a prescribed composition). Fig. [3] is the probability 
distribution of the volume at coexistence, the two equal height peaks of the distribution 
correspond to the coexisting fluid and solid phases, respectively. In Fig HI we plot the phase 
diagram of the polydisperse HSS in the {p,S) plane. For the polydispersity used in the 
simulation, the phase diagram has the same properties as the previous simulation 3|. The 
transition points reported here are not accurately located because extensive runs are needed 
to refine the locations, but it is suffice to exhibit the validity of our method. 

To conclude we provide a new Monte Carlo method to study the fluid-solid phase tran- 
sition of the polydisperse colloidal system. The validity of the method is demonstrated by 
locating the fluid-solid transition point of the polydisperse HHS. The obtained results are 
consistent with the previous results. The method can be directly extended to the poly- 
disperse soft-interaction systems and the hard-interaction system with other polydisperse 
attributes. It presents a complement to the current Gibbs-Duhem integration method [Sj. 
The estimate of cloud points is important to fully understand the phase behavior of polydis- 
perse colloids [30|. It is possible to extend our current method to determine the cloud points 
for the fluid-solid transition of a polydisperse system. We are working on this direction and 
the results will be reported in the near future. 
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